## Replication code for: Katsumata, Hiroto and Shunya Noda. 2023. 
## "'Kick Them Out' as a Voting Strategy: Theory and Evidence from Multi-member District Elections." 
## The Journal of Politics, forthcoming.

## Re-analysis of Hix, Hortala-Vallve, and Riambau-Armet (2016) JoP.
## Author: Hiroto Katsumata
## Date: May, 2023

## Double check that you have not load tidyverse package yet
detach(package:tidyverse,unload=TRUE)
# You will get an error message but it is okay

## Load packages
library(Hmisc)
library(MASS)
library(sandwich)
library(tidyverse)

## Load functions
source("functions/kickout_functions.R")


## Data wrangling ====================

## Read and wrangle data
data <- read.csv("data/Data_DistrictMagnitudeAugust2016.csv", stringsAsFactors = FALSE) %>% 
				filter(M %in% 2:3) %>%
				mutate(p_cluster = ceiling(period / 5),
							 cluster = subject * 1000 + session * 100 + p_cluster)

## M == 2
dataM2 <- data %>% filter(M == 2) %>%
					mutate(cluster = factor(cluster))
dataM2 <- dataM2 %>% 
					dplyr::select(L_VOTE_1, L_VOTE_2, L_VOTE_3, L_VOTE_4, L_VOTE_5) %>%
					as.matrix() %>%
					saintelague(lvote = ., m = 2) %>%
					data.frame(dataM2, .) %>%
					filter(n_topmplus1 == 3)
dataM2 <- preftop(data = dataM2) %>%
					data.frame(dataM2, .) %>%
					mutate(vote_runnerup = vote == runnerup,
								 vote_topmplus1_first = vote == topmplus1_first,
								 vote_topmplus1_second = vote == topmplus1_second,
								 vote_notrunuptopm_first = vote == notrunuptopm_first,
								 vote_notrunuptopm_second = vote == notrunuptopm_second,
								 vote_lastwin = vote == ranksecond) %>%
					mutate(pref_notrunuptopm_least = pref_notrunuptopm_second) %>%
					filter(nondiscriminate == 0)

dataM2$vote_notrunuptopm <- NA
for (i in 1:nrow(dataM2)) {
	dataM2$vote_notrunuptopm[i] <- dataM2$vote[i] %in% 
														c(1:5)[as.logical(dataM2[i, c("notrunuptopm_1", "notrunuptopm_2",
								 												"notrunuptopm_3", "notrunuptopm_4", "notrunuptopm_5")])]
}

## M == 3
dataM3 <- data %>% filter(M == 3) %>%
					mutate(cluster = factor(cluster))
dataM3 <- dataM3 %>% 
					select(L_VOTE_1, L_VOTE_2, L_VOTE_3, L_VOTE_4, L_VOTE_5) %>%
					as.matrix() %>%
					saintelague(lvote = ., m = 3) %>%
					data.frame(dataM3, .)
dataM3 <- preftop(data = dataM3) %>%
					data.frame(dataM3, .) %>%
					mutate(vote_runnerup = vote == runnerup,
								 vote_topmplus1_first = vote == topmplus1_first,
								 vote_topmplus1_second = vote == topmplus1_second,
								 vote_notrunuptopm_first = vote == notrunuptopm_first,
								 vote_notrunuptopm_second = vote == notrunuptopm_second,
								 vote_lastwin = vote == rankthird) %>%
					mutate(pref_notrunuptopm_least = ifelse(n_topmplus1 == 4, pref_notrunuptopm_third, 
																																		pref_notrunuptopm_second),
								 pref_topmplus1_least = ifelse(n_topmplus1 == 4, pref_topmplus1_fourth, 
																																 pref_topmplus1_third)) %>%
					filter(nondiscriminate == 0)

dataM3$vote_notrunuptopm <- NA
for (i in 1:nrow(dataM3)) {
	dataM3$vote_notrunuptopm[i] <- dataM3$vote[i] %in% 
														c(1:5)[as.logical(dataM3[i, c("notrunuptopm_1", "notrunuptopm_2",
								 												"notrunuptopm_3", "notrunuptopm_4", "notrunuptopm_5")])]
}


## Select variables used in the analysis
dataM2 <- dataM2 %>%
					dplyr::select(c(p_cluster, cluster, n_topmplus1, 
													pref_topmplus1_first, pref_topmplus1_second, pref_topmplus1_third, pref_topmplus1_fourth, 
													pref_runnerup,
													pref_notrunuptopm_first, pref_notrunuptopm_second, pref_notrunuptopm_third, pref_notrunuptopm_least, 
													vote_runnerup, vote_lastwin,
													vote_topmplus1_first, vote_topmplus1_second, 
													vote_notrunuptopm_first, vote_notrunuptopm_second, vote_notrunuptopm))

dataM3 <- dataM3 %>%
					dplyr::select(c(p_cluster, cluster, n_topmplus1, 
													pref_topmplus1_first, pref_topmplus1_second, pref_topmplus1_third, pref_topmplus1_fourth, 
													pref_topmplus1_least,
													pref_runnerup, 
													pref_notrunuptopm_first, pref_notrunuptopm_second, pref_notrunuptopm_third, pref_notrunuptopm_least, 
													vote_runnerup, vote_lastwin,
													vote_topmplus1_first, vote_topmplus1_second, 
													vote_notrunuptopm_first, vote_notrunuptopm_second, vote_notrunuptopm))


## Estimation ====================

## Main data for the analysis: "dataM2" and "dataM3"
## See the codebook "kickout_codebook.txt" for details
head(dataM2)
head(dataM3)


## A: Kick out the least preferred party among the top M + 1 parties
## Vote for the second-most preferred party against the most preferred one
##  among the top M + 1 parties
## Treatment: preference for the least preferred party among the viable parties

## M == 2
## All Periods
result_M2A <- 
subset(dataM2, vote_topmplus1_first + vote_topmplus1_second == 1) %>%
summary_cluster_a(data = ., formula = vote_topmplus1_second ~ 
																			pref_topmplus1_first + 
																			pref_topmplus1_second + 
																			pref_topmplus1_third,
														type = "a")

## Periods 1--6
result_M2A_1 <- 
subset(dataM2, vote_topmplus1_first + vote_topmplus1_second == 1) %>% 
filter(p_cluster <= 6) %>%
summary_cluster_a(data = ., formula = vote_topmplus1_second ~ 
																			pref_topmplus1_first + 
																			pref_topmplus1_second + 
																			pref_topmplus1_third,
														type = "a")
## Periods 4--9
result_M2A_2 <- 
subset(dataM2, vote_topmplus1_first + vote_topmplus1_second == 1) %>% 
filter(p_cluster >= 4 & p_cluster <= 9) %>%
summary_cluster_a(data = ., formula = vote_topmplus1_second ~ 
																			pref_topmplus1_first + 
																			pref_topmplus1_second + 
																			pref_topmplus1_third,
														type = "a")
## Periods 7--12
result_M2A_3 <- 
subset(dataM2, vote_topmplus1_first + vote_topmplus1_second == 1) %>% 
filter(p_cluster >= 7 & p_cluster <= 12) %>%
summary_cluster_a(data = ., formula = vote_topmplus1_second ~ 
																			pref_topmplus1_first + 
																			pref_topmplus1_second + 
																			pref_topmplus1_third,
														type = "a")

## M == 3
## All Periods
result_M3A <- 
subset(dataM3, vote_topmplus1_first + vote_topmplus1_second == 1) %>%
summary_cluster_a(data = ., formula = vote_topmplus1_second ~ 
																			pref_topmplus1_first + 
																			pref_topmplus1_second + 
																			pref_topmplus1_least +
																			n_topmplus1,
														type = "a")
## Periods 1--6
result_M3A_1 <- 
subset(dataM3, vote_topmplus1_first + vote_topmplus1_second == 1) %>% 
filter(p_cluster <= 6) %>%
summary_cluster_a(data = ., formula = vote_topmplus1_second ~ 
																			pref_topmplus1_first + 
																			pref_topmplus1_second + 
																			pref_topmplus1_least +
																			n_topmplus1,
														type = "a")
## Periods 4--9
result_M3A_2 <- 
subset(dataM3, vote_topmplus1_first + vote_topmplus1_second == 1) %>% 
filter(p_cluster >= 4 & p_cluster <= 9) %>%
summary_cluster_a(data = ., formula = vote_topmplus1_second ~ 
																			pref_topmplus1_first + 
																			pref_topmplus1_second + 
																			pref_topmplus1_least +
																			n_topmplus1,
														type = "a")
## Periods 7--12
result_M3A_3 <- 
subset(dataM3, vote_topmplus1_first + vote_topmplus1_second == 1) %>% 
filter(p_cluster >= 7 & p_cluster <= 12) %>%
summary_cluster_a(data = ., formula = vote_topmplus1_second ~ 
																			pref_topmplus1_first + 
																			pref_topmplus1_second + 
																			pref_topmplus1_least +
																			n_topmplus1,
														type = "a")

## B: Kick out the runner-up
## Vote for the second-most preferred party against the most preferred one
##  among the top M parties excluding the runner-up party
## Treatment: preference for the runner-up

## M == 2
## All Periods
result_M2B <- 
subset(dataM2, vote_notrunuptopm_first + vote_notrunuptopm_second == 1) %>%
summary_cluster_a(data = ., formula = vote_notrunuptopm_second ~ 
																			pref_notrunuptopm_first + 
																			pref_notrunuptopm_second + 
																			pref_runnerup,
														type = "a")
## Periods 1--6
result_M2B_1 <- 
subset(dataM2, vote_notrunuptopm_first + vote_notrunuptopm_second == 1) %>%
filter(p_cluster <= 6) %>%
summary_cluster_a(data = ., formula = vote_notrunuptopm_second ~ 
																			pref_notrunuptopm_first + 
																			pref_notrunuptopm_second + 
																			pref_runnerup,
														type = "a")
## Periods 4--9
result_M2B_2 <- 
subset(dataM2, vote_notrunuptopm_first + vote_notrunuptopm_second == 1) %>%
filter(p_cluster >= 4 & p_cluster <= 9) %>%
summary_cluster_a(data = ., formula = vote_notrunuptopm_second ~ 
																			pref_notrunuptopm_first + 
																			pref_notrunuptopm_second + 
																			pref_runnerup,
														type = "a")
## Periods 7--12
result_M2B_3 <- 
subset(dataM2, vote_notrunuptopm_first + vote_notrunuptopm_second == 1) %>%
filter(p_cluster >= 7 & p_cluster <= 12) %>%
summary_cluster_a(data = ., formula = vote_notrunuptopm_second ~ 
																			pref_notrunuptopm_first + 
																			pref_notrunuptopm_second + 
																			pref_runnerup,
														type = "a")

## M == 3
## All Periods
result_M3B <- 
subset(dataM3, vote_notrunuptopm_first + vote_notrunuptopm_second == 1) %>%
summary_cluster_a(data = ., formula = vote_notrunuptopm_second ~ 
																			pref_notrunuptopm_first + 
																			pref_notrunuptopm_second + 
																			pref_runnerup +
																			n_topmplus1,
														type = "a")
## Periods 1--6
result_M3B_1 <- 
subset(dataM3, vote_notrunuptopm_first + vote_notrunuptopm_second == 1) %>%
filter(p_cluster <= 6) %>%
summary_cluster_a(data = ., formula = vote_notrunuptopm_second ~ 
																			pref_notrunuptopm_first + 
																			pref_notrunuptopm_second + 
																			pref_runnerup +
																			n_topmplus1,
														type = "a")
## Periods 4--9
result_M3B_2 <- 
subset(dataM3, vote_notrunuptopm_first + vote_notrunuptopm_second == 1) %>%
filter(p_cluster >= 4 & p_cluster <= 9) %>%
summary_cluster_a(data = ., formula = vote_notrunuptopm_second ~ 
																			pref_notrunuptopm_first + 
																			pref_notrunuptopm_second + 
																			pref_runnerup +
																			n_topmplus1,
														type = "a")
## Periods 7--12
result_M3B_3 <- 
subset(dataM3, vote_notrunuptopm_first + vote_notrunuptopm_second == 1) %>%
filter(p_cluster >= 7 & p_cluster <= 12) %>%
summary_cluster_a(data = ., formula = vote_notrunuptopm_second ~ 
																			pref_notrunuptopm_first + 
																			pref_notrunuptopm_second + 
																			pref_runnerup +
																			n_topmplus1,
														type = "a")

## C: Kick out the least preferred party among the top M parties
## Vote for the runner-up party against the top M parties
## Treatment: preference for the least preferred party among the top M parties

## M == 2
## All Periods
result_M2C <- 
subset(dataM2, (vote_notrunuptopm == 1 | vote_runnerup == 1)) %>%
summary_cluster_a(data = ., formula = vote_runnerup ~ 
																			pref_runnerup + pref_notrunuptopm_first + 
																			pref_notrunuptopm_least,
														type = "a")
## Exclude voters who vote for the Mth-place party
result_M2Cr <- 
subset(dataM2, (vote_notrunuptopm == 1 | vote_runnerup == 1) & vote_lastwin == 0) %>%
summary_cluster_a(data = ., formula = vote_runnerup ~ 
																			pref_runnerup + pref_notrunuptopm_first + 
																			pref_notrunuptopm_least,
														type = "a")
## Periods 1--6
result_M2C_1 <- 
subset(dataM2, (vote_notrunuptopm == 1 | vote_runnerup == 1)) %>%
filter(p_cluster <= 6) %>%
summary_cluster_a(data = ., formula = vote_runnerup ~ 
																			pref_runnerup + pref_notrunuptopm_first + 
																			pref_notrunuptopm_least,
														type = "a")
## Periods 4--9
result_M2C_2 <- 
subset(dataM2, (vote_notrunuptopm == 1 | vote_runnerup == 1)) %>%
filter(p_cluster >= 4 & p_cluster <= 9) %>%
summary_cluster_a(data = ., formula = vote_runnerup ~ 
																			pref_runnerup + pref_notrunuptopm_first + 
																			pref_notrunuptopm_least,
														type = "a")
## Periods 7--12
result_M2C_3 <- 
subset(dataM2, (vote_notrunuptopm == 1 | vote_runnerup == 1)) %>%
filter(p_cluster >= 7 & p_cluster <= 12) %>%
summary_cluster_a(data = ., formula = vote_runnerup ~ 
																			pref_runnerup + pref_notrunuptopm_first + 
																			pref_notrunuptopm_least,
														type = "a")

## M == 3
## All Periods
result_M3C <- 
subset(dataM3, (vote_notrunuptopm == 1 | vote_runnerup == 1)) %>%
summary_cluster_a(data = ., formula = vote_runnerup ~ 
																			pref_runnerup + pref_notrunuptopm_first + 
																			pref_notrunuptopm_second + pref_notrunuptopm_least +
																			n_topmplus1,
														type = "b2")
## Exclude voters who vote for the Mth-place party
result_M3Cr <- 
subset(dataM3, (vote_notrunuptopm == 1 | vote_runnerup == 1) & vote_lastwin == 0) %>%
summary_cluster_a(data = ., formula = vote_runnerup ~ 
																			pref_runnerup + pref_notrunuptopm_first + 
																			pref_notrunuptopm_second + pref_notrunuptopm_least +
																			n_topmplus1,
														type = "b2")
## Periods 1--6
result_M3C_1 <- 
subset(dataM3, (vote_notrunuptopm == 1 | vote_runnerup == 1)) %>%
filter(p_cluster <= 6) %>%
summary_cluster_a(data = ., formula = vote_runnerup ~ 
																			pref_runnerup + pref_notrunuptopm_first + 
																			pref_notrunuptopm_second + pref_notrunuptopm_least +
																			n_topmplus1,
														type = "b2")
## Periods 4--9
result_M3C_2 <- 
subset(dataM3, (vote_notrunuptopm == 1 | vote_runnerup == 1)) %>%
filter(p_cluster >= 4 & p_cluster <= 9) %>%
summary_cluster_a(data = ., formula = vote_runnerup ~ 
																			pref_runnerup + pref_notrunuptopm_first + 
																			pref_notrunuptopm_second + pref_notrunuptopm_least +
																			n_topmplus1,
														type = "b2")
## Periods 7--12
result_M3C_3 <- 
subset(dataM3, (vote_notrunuptopm == 1 | vote_runnerup == 1)) %>%
filter(p_cluster >= 7 & p_cluster <= 12) %>%
summary_cluster_a(data = ., formula = vote_runnerup ~ 
																			pref_runnerup + pref_notrunuptopm_first + 
																			pref_notrunuptopm_second + pref_notrunuptopm_least +
																			n_topmplus1,
														type = "b2")

## Visualize the average marginal effects
## A: Kick out the least preferred party among the top M + 1 parties
## M == 3
xlabel <- "Preference for the least preferred viable party"
ylabel <- "Prob. of voting for the second-most preferred viable party        \nvs. most preferred viable party"
subset(dataM3, vote_topmplus1_first + vote_topmplus1_second == 1)$pref_topmplus1_least %>%
summary()
set.seed(20201234)
sim_data <- glm(data = subset(dataM3, 
															vote_topmplus1_first + vote_topmplus1_second == 1), 
								formula = vote_topmplus1_second ~ 
													pref_topmplus1_least + 
													pref_topmplus1_first + 
													pref_topmplus1_second +
													n_topmplus1, 
								family = binomial(probit),
								x = TRUE)$x
res <- glm_cluster(data = subset(dataM3, 
																 vote_topmplus1_first + vote_topmplus1_second == 1), 
									 formula = vote_topmplus1_second ~ 
														 pref_topmplus1_least + 
														 pref_topmplus1_first + 
														 pref_topmplus1_second +
														 n_topmplus1, 
									 cluster = "cluster", 
									 family = binomial(probit)) %>%
			 ame_cluster(object = ., 
			 		 				 sim_data = sim_data, 
			 						 treat = seq(0, 150, by = 5), B = 10000, 
			 						 xlab = xlabel,
			 						 ylab = ylabel,
			 						 ylab2 = "Number of observations",
			 						 y2scale = 3,
			 						 histobreaks = seq(-2.5, 152.5, by = 5))
## Figure C.1
res$figure
ggsave("figures/experimentM3A_ame.pdf", width = 19, height = 19, units = "cm")
q13 <- subset(dataM3, vote_topmplus1_first + vote_topmplus1_second == 1)$pref_topmplus1_least %>%
			 quantile(., probs = c(0.25, 0.75))
res <- glm_cluster(data = subset(dataM3, 
																 vote_topmplus1_first + vote_topmplus1_second == 1), 
									 formula = vote_topmplus1_second ~ 
														 pref_topmplus1_least + 
														 pref_topmplus1_first + 
														 pref_topmplus1_second +
														 n_topmplus1, 
									 cluster = "cluster", 
									 family = binomial(probit)) %>%
			 amedif_cluster(object = ., 
			 		 						sim_data = sim_data, 
			 								treat = c(q13[1], q13[2]), B = 10000)
res

## B: Kick out the runner-up
## M == 3
xlabel <- "Preference for the runner-up party"
ylabel <- "Prob. of voting for the second-most preferred    \n leading party vs. most preferred leading party"
subset(dataM3, vote_notrunuptopm_first + vote_notrunuptopm_second == 1)$pref_runnerup %>%
summary()
set.seed(20201234)
sim_data <- glm(data = subset(dataM3, vote_notrunuptopm_first + vote_notrunuptopm_second == 1), 
								formula = vote_notrunuptopm_second ~ 
													pref_runnerup +
													pref_notrunuptopm_first + 
													pref_notrunuptopm_second +
													n_topmplus1, 
								family = binomial(probit),
								x = TRUE)$x
res <- glm_cluster(data = subset(dataM3, vote_notrunuptopm_first + vote_notrunuptopm_second == 1),
									 formula = vote_notrunuptopm_second ~ 
														 pref_runnerup +
														 pref_notrunuptopm_first + 
														 pref_notrunuptopm_second +
														n_topmplus1, 
									 cluster = "cluster", 
									 family = binomial(probit)) %>%
			 ame_cluster(object = ., 
			 		 				 sim_data = sim_data, 
			 						 treat = seq(0, 190, by = 5), B = 10000, 
			 						 xlab = xlabel,
			 						 ylab = ylabel,
			 						 ylab2 = "Number of observations",
			 						 y2scale = 3,
			 						 histobreaks = seq(-2.5, 192.5, by = 5))
## Figure 4
res$figure
ggsave("figures/experimentM3B_ame.pdf", width = 19, height = 19, units = "cm")
q13 <- subset(dataM3, vote_notrunuptopm_first + vote_notrunuptopm_second == 1)$pref_runnerup %>%
			 quantile(., probs = c(0.25, 0.75))
res <- glm_cluster(data = subset(dataM3, vote_notrunuptopm_first + vote_notrunuptopm_second == 1), 
									 formula = vote_notrunuptopm_second ~ 
														 pref_runnerup +
														 pref_notrunuptopm_first + 
														 pref_notrunuptopm_second +
														 n_topmplus1, 
									 cluster = "cluster", 
									 family = binomial(probit)) %>%
			 amedif_cluster(object = ., 
			 		 						sim_data = sim_data, 
			 								treat = c(q13[1], q13[2]), B = 10000)
res

## Comparison to the effect of the preference for the most preferred leading party
set.seed(20201234)
sim_data_compare <- glm(data = subset(dataM3, vote_notrunuptopm_first + vote_notrunuptopm_second == 1), 
								formula = vote_notrunuptopm_second ~ 
													pref_notrunuptopm_first + 
													pref_runnerup +
													pref_notrunuptopm_second +
													n_topmplus1, 
								family = binomial(probit),
								x = TRUE)$x
q13_compare <- subset(dataM3, vote_notrunuptopm_first + vote_notrunuptopm_second == 1)$pref_notrunuptopm_first %>%
			 quantile(., probs = c(0.25, 0.75))
res <- glm_cluster(data = subset(dataM3, vote_notrunuptopm_first + vote_notrunuptopm_second == 1), 
									 formula = vote_notrunuptopm_second ~ 
														 pref_notrunuptopm_first + 
														 pref_runnerup +
														 pref_notrunuptopm_second +
														 n_topmplus1, 
									 cluster = "cluster", 
									 family = binomial(probit)) %>%
			 amedif_cluster(object = ., 
			 		 						sim_data = sim_data_compare, 
			 								treat = c(q13_compare[1], q13_compare[2]), B = 10000)
res

## Comparison to the effect of the preference for the second-most preferred leading party
set.seed(20201234)
sim_data_compare <- glm(data = subset(dataM3, vote_notrunuptopm_first + vote_notrunuptopm_second == 1), 
								formula = vote_notrunuptopm_second ~ 
													pref_notrunuptopm_second + 
													pref_runnerup +
													pref_notrunuptopm_first +
													n_topmplus1, 
								family = binomial(probit),
								x = TRUE)$x
q13_compare <- subset(dataM3, vote_notrunuptopm_first + vote_notrunuptopm_second == 1)$pref_notrunuptopm_second %>%
			 quantile(., probs = c(0.25, 0.75))
res <- glm_cluster(data = subset(dataM3, vote_notrunuptopm_first + vote_notrunuptopm_second == 1), 
									 formula = vote_notrunuptopm_second ~ 
														 pref_notrunuptopm_second + 
														 pref_runnerup +
														 pref_notrunuptopm_first +
														 n_topmplus1, 
									 cluster = "cluster", 
									 family = binomial(probit)) %>%
			 amedif_cluster(object = ., 
			 		 						sim_data = sim_data_compare, 
			 								treat = c(q13_compare[1], q13_compare[2]), B = 10000)
res


## C: Kick out the least preferred party within the top M parties
## M == 2
xlabel <- "    Preference for the second-most preferred leading party"
ylabel <- "Prob. of voting for the runner-up vs. the leading parties    "
subset(dataM2, (vote_notrunuptopm == 1 | vote_runnerup == 1))$pref_notrunuptopm_least %>%
summary()
set.seed(20201234)
sim_data <- glm(data = subset(dataM2, 
															(vote_notrunuptopm == 1 | vote_runnerup == 1)), 
								formula = vote_runnerup ~ 
													pref_notrunuptopm_least +
													pref_runnerup + 
													pref_notrunuptopm_first, 
								family = binomial(probit),
								x = TRUE)$x
res <- glm_cluster(data = subset(dataM2, 
																 (vote_notrunuptopm == 1 | vote_runnerup == 1)), 
									 formula = vote_runnerup ~ 
														 pref_notrunuptopm_least +
														 pref_runnerup + 
														 pref_notrunuptopm_first, 
									 cluster = "cluster", 
									 family = binomial(probit)) %>%
			 ame_cluster(object = ., 
			 		 				 sim_data = sim_data, 
			 						 treat = seq(30, 230, by = 5), B = 10000, 
			 						 xlab = xlabel,
			 						 ylab = ylabel,
			 						 ylab2 = "Number of observations",
			 						 y2scale = 3,
			 						 histobreaks = seq(27.5, 232.5, by = 5))
res$figure
## Figure 3
ggsave("figures/experimentM2C_ame.pdf", width = 19, height = 19, units = "cm")
q13 <- subset(dataM2, (vote_notrunuptopm == 1 | vote_runnerup == 1))$pref_notrunuptopm_least %>%
			 quantile(., probs = c(0.25, 0.75))
res <- glm_cluster(data = subset(dataM2, 
																 (vote_notrunuptopm == 1 | vote_runnerup == 1)), 
									 formula = vote_runnerup ~ 
														 pref_notrunuptopm_least +
														 pref_runnerup + 
														 pref_notrunuptopm_first, 
									 cluster = "cluster", 
									 family = binomial(probit)) %>%
			 amedif_cluster(object = ., 
			 		 				 		sim_data = sim_data, 
			 								treat = c(q13[1], q13[2]), B = 10000)
res

## Comparison to the effect of the preference for the runner-up party
set.seed(20201234)
sim_data_compare <- glm(data = subset(dataM2, 
															(vote_notrunuptopm == 1 | vote_runnerup == 1)), 
								formula = vote_runnerup ~ 
													pref_runnerup + 
													pref_notrunuptopm_least +
													pref_notrunuptopm_first, 
								family = binomial(probit),
								x = TRUE)$x
q13_compare <- subset(dataM2, (vote_notrunuptopm == 1 | vote_runnerup == 1))$pref_runnerup %>%
			 quantile(., probs = c(0.25, 0.75))
res <- glm_cluster(data = subset(dataM2, 
																 (vote_notrunuptopm == 1 | vote_runnerup == 1)), 
									 formula = vote_runnerup ~ 
														 pref_runnerup + 
														 pref_notrunuptopm_least +
														 pref_notrunuptopm_first, 
									 cluster = "cluster", 
									 family = binomial(probit)) %>%
			 amedif_cluster(object = ., 
			 		 				 		sim_data = sim_data_compare, 
			 								treat = c(q13_compare[1], q13_compare[2]), B = 10000)
res



## Tables
## LaTeX table
## A: Kick out the least preferred party among the top M + 1 parties
## Table C.1
rnameA <- c("District magnitude ($M$)",
						"Pref. for the most preferred",
						"~~viable party",
						"Pref. for the second-most preferred",
						"~~viable party",
						"Pref. for the least preferred",
						"~~viable party",
						"Number of observations",
						"Number of clusters")
resultA <- data.frame(rnameA,
											NA,
											c("2", result_M2A),
											NA,
											c("3", result_M3A))
clabelsA <- c("", "", "Model A1", "", "Model A2")
captionA <- "(Lab Experiment) Kicking out the Least Preferred Viable Party"
noteA <- 
"\\parbox{0.7\\textwidth}
{Notes: Coefficients and standard errors are multiplied by 100. 
The outcome variable is voting for the second-most preferred viable party
against most preferred viable party.
The model for $M = 3$ includes the number of viable parties as covariates.
Cluster-robust standard errors in parentheses.}"
latex(object = resultA,
      title = "",
      file = "tables/experimentA.tex",
      label = "tb_experimentA",
      caption = captionA,
      insert.bottom = noteA,
      first.hline.double = FALSE,
      rowname = NULL,
      colheads = clabelsA,
      col.just = c("l", rep("c", ncol(resultA) - 1)),
      n.rgroup = c(1, 6, 2),
      longtable = FALSE,
      center = "centering")

## Learning
## M = 2
resultA <- data.frame(rnameA,
											NA,
											c("2", result_M2A_1),
											NA,
											c("2", result_M2A_2),
											NA,
											c("2", result_M2A_3),
											NA,
											c("2", result_M2A))
clabelsPeriod <- c("", "", "Period 1--6", "", "Period 4--9", "", "Period 7--12", "", "All Periods")
captionA_M2 <- "(Lab Experiment) Learning Effect: Kicking out the Least Preferred Viable Party ($M = 2$)"
noteA_M2 <- 
"\\parbox{0.98\\textwidth}
{Notes: Coefficients and standard errors are multiplied by 100. 
The outcome variable is voting for the second-most preferred viable party
against most preferred viable party.
Cluster-robust standard errors in parentheses.}"
latex(object = resultA,
      title = "",
      file = "tables/experimentA_learning_M2.tex",
      label = "tb_experimentA_learning_M2",
      caption = captionA_M2,
      insert.bottom = noteA_M2,
      first.hline.double = FALSE,
      rowname = NULL,
      colheads = clabelsPeriod,
      col.just = c("l", rep("c", ncol(resultA) - 1)),
      n.rgroup = c(1, 6, 2),
      longtable = FALSE,
      center = "centering")

## M = 3
resultA <- data.frame(rnameA,
											NA,
											c("3", result_M3A_1),
											NA,
											c("3", result_M3A_2),
											NA,
											c("3", result_M3A_3),
											NA,
											c("3", result_M3A))
captionA_M3 <- "(Lab Experiment) Learning Effect: Kicking out the Least Preferred Viable Party ($M = 3$)"
noteA_M3 <- 
"\\parbox{0.98\\textwidth}
{Notes: Coefficients and standard errors are multiplied by 100. 
The outcome variable is voting for the second-most preferred viable party
against most preferred viable party.
The number of viable parties is included as covariates.
Cluster-robust standard errors in parentheses.}"
latex(object = resultA,
      title = "",
      file = "tables/experimentA_learning_M3.tex",
      label = "tb_experimentA_learning_M3",
      caption = captionA_M3,
      insert.bottom = noteA_M3,
      first.hline.double = FALSE,
      rowname = NULL,
      colheads = clabelsPeriod,
      col.just = c("l", rep("c", ncol(resultA) - 1)),
      n.rgroup = c(1, 6, 2),
      longtable = FALSE,
      center = "centering")


## B: Kick out the runner-up
## Table D.3
rnameB <- c("District magnitude ($M$)",
						"Pref. for the most preferred",
						"~~leading party",
						"Pref. for the second-most preferred",
						"~~leading party",
						"\\multirow{2}{*}{Pref. for the runner-up party}",
						"",
						"Number of observations",
						"Number of clusters")
resultB <- data.frame(rnameB,
											NA,
											c("2", result_M2B),
											NA,
											c("3", result_M3B))
clabelsB <- c("", "", "Model B1", "", "Model B2")
captionB <- "(Lab Experiment) Kicking out the Runner-up Party"
noteB <- 
"\\parbox{0.7\\textwidth}
{Notes: Coefficients and standard errors are multiplied by 100. 
The outcome variable is voting for the second-most preferred leading party
against most preferred leading party.
The model for $M = 3$ includes the number of viable parties as covariates.
Cluster-robust standard errors in parentheses.}"
latex(object = resultB,
      title = "",
      file = "tables/experimentB.tex",
      label = "tb_experimentB",
      caption = captionB,
      insert.bottom = noteB,
      first.hline.double = FALSE,
      rowname = NULL,
      colheads = clabelsB,
      col.just = c("l", rep("c", ncol(resultB) - 1)),
      n.rgroup = c(1, 6, 2),
      longtable = FALSE,
      center = "centering")

## Learning
## M = 2
## Table E.3
resultB <- data.frame(rnameB,
											NA,
											c("2", result_M2B_1),
											NA,
											c("2", result_M2B_2),
											NA,
											c("2", result_M2B_3),
											NA,
											c("2", result_M2B))
captionB_M2 <- "(Lab Experiment) Learning Effect: Kicking out the Runner-up Party ($M = 2$)"
noteB_M2 <- 
"\\parbox{0.98\\textwidth}
{Notes: Coefficients and standard errors are multiplied by 100. 
The outcome variable is voting for the second-most preferred leading party
against most preferred leading party.
Cluster-robust standard errors in parentheses.}"
latex(object = resultB,
      title = "",
      file = "tables/experimentB_learning_M2.tex",
      label = "tb_experimentB_learning_M2",
      caption = captionB_M2,
      insert.bottom = noteB_M2,
      first.hline.double = FALSE,
      rowname = NULL,
      colheads = clabelsPeriod,
      col.just = c("l", rep("c", ncol(resultB) - 1)),
      n.rgroup = c(1, 6, 2),
      longtable = FALSE,
      center = "centering")

## M = 3
## Table E.4
resultB <- data.frame(rnameB,
											NA,
											c("3", result_M3B_1),
											NA,
											c("3", result_M3B_2),
											NA,
											c("3", result_M3B_3),
											NA,
											c("3", result_M3B))
captionB_M3 <- "(Lab Experiment) Learning Effect: Kicking out the Runner-up Party ($M = 3$)"
noteB_M3 <- 
"\\parbox{0.98\\textwidth}
{Notes: Coefficients and standard errors are multiplied by 100. 
The outcome variable is voting for the second-most preferred leading party
against most preferred leading party.
The number of viable parties is included as covariates.
Cluster-robust standard errors in parentheses.}"
latex(object = resultB,
      title = "",
      file = "tables/experimentB_learning_M3.tex",
      label = "tb_experimentB_learning_M3",
      caption = captionB_M3,
      insert.bottom = noteB_M3,
      first.hline.double = FALSE,
      rowname = NULL,
      colheads = clabelsPeriod,
      col.just = c("l", rep("c", ncol(resultB) - 1)),
      n.rgroup = c(1, 6, 2),
      longtable = FALSE,
      center = "centering")


## C: Kick out the least preferred party within top M parties
## Table D.2
rnameC <- c("District magnitude ($M$)",
						"\\multirow{2}{*}{Pref. for the runner-up party}",
						"",
						"Pref. for the most preferred",
						"~~leading party",
						"Pref. for the second-most preferred",
						"~~leading party",
						"Pref. for the third-most preferred",
						"~~leading party",
						"Number of observations",
						"Number of clusters")
resultC <- data.frame(rnameC,
											NA,
											c("2", result_M2C[1:6], NA, NA, result_M2C[7:8]),
											NA,
											c("3", result_M3C),
											NA,
											c("2", result_M2Cr[1:6], NA, NA, result_M2Cr[7:8]),
											NA,
											c("3", result_M3Cr))
clabelsC <- c("", "", "Model C1", "", "Model C2", "", "Model C3", "", "Model C4")
captionC <- "(Lab Experiment) Kicking out the Least Preferred Leading Party"
noteC <- 
"\\parbox{0.98\\textwidth}
{Notes: Coefficients and standard errors are multiplied by 100. 
The outcome variable is voting for the runner-up
against one of the leading parties.
The models for $M = 3$ include the number of viable parties as covariates.
The samples for Models C3 and C4 exclude voters who vote for the $M$th-place party.
Cluster-robust standard errors in parentheses.}"
latex(object = resultC,
      title = "",
      file = "tables/experimentC.tex",
      label = "tb_experimentC",
      caption = captionC,
      insert.bottom = noteC,
      first.hline.double = FALSE,
      rowname = NULL,
      colheads = clabelsC,
      col.just = c("l", rep("c", ncol(resultC) - 1)),
      n.rgroup = c(1, 8, 2),
      longtable = FALSE,
      center = "centering")

## Learning
## M = 2
## Table E.1
rnameC_M2 <- c("District magnitude ($M$)",
						"\\multirow{2}{*}{Pref. for the runner-up party}",
						"",
						"Pref. for the most preferred",
						"~~leading party",
						"Pref. for the second-most preferred",
						"~~leading party",
						"Number of observations",
						"Number of clusters")
resultC <- data.frame(rnameC_M2,
											NA,
											c("2", result_M2C_1),
											NA,
											c("2", result_M2C_2),
											NA,
											c("2", result_M2C_3),
											NA,
											c("2", result_M2C))
captionC_M2 <- "(Lab Experiment) Learning Effect: Kicking out the Least Preferred Leading Party ($M = 2$)"
noteC_M2 <- 
"\\parbox{0.98\\textwidth}
{Notes: Coefficients and standard errors are multiplied by 100. 
The outcome variable is voting for the runner-up
against one of the leading parties.
Cluster-robust standard errors in parentheses.}"
latex(object = resultC,
      title = "",
      file = "tables/experimentC_learning_M2.tex",
      label = "tb_experimentC_learning_M2",
      caption = captionC_M2,
      insert.bottom = noteC_M2,
      first.hline.double = FALSE,
      rowname = NULL,
      colheads = clabelsPeriod,
      col.just = c("l", rep("c", ncol(resultC) - 1)),
      n.rgroup = c(1, 6, 2),
      longtable = FALSE,
      center = "centering")

## M = 3
## Table E.2
rnameC_M3 <- c("District magnitude ($M$)",
						"\\multirow{2}{*}{Pref. for the runner-up party}",
						"",
						"Pref. for the most preferred",
						"~~leading party",
						"Pref. for the second-most preferred",
						"~~leading party",
						"Pref. for the third-most preferred",
						"~~leading party",
						"Number of observations",
						"Number of clusters")
resultC <- data.frame(rnameC_M3,
											NA,
											c("3", result_M3C_1),
											NA,
											c("3", result_M3C_2),
											NA,
											c("3", result_M3C_3),
											NA,
											c("3", result_M3C))
captionC_M3 <- "(Lab Experiment) Learning Effect: Kicking out the Least Preferred Leading Party ($M = 3$)"
noteC_M3 <- 
"\\parbox{0.98\\textwidth}
{Notes: Coefficients and standard errors are multiplied by 100. 
The outcome variable is voting for the runner-up
against one of the leading parties.
The number of viable parties is included as covariates.
Cluster-robust standard errors in parentheses.}"
latex(object = resultC,
      title = "",
      file = "tables/experimentC_learning_M3.tex",
      label = "tb_experimentC_learning_M3",
      caption = captionC_M3,
      insert.bottom = noteC_M3,
      first.hline.double = FALSE,
      rowname = NULL,
      colheads = clabelsPeriod,
      col.just = c("l", rep("c", ncol(resultC) - 1)),
      n.rgroup = c(1, 8, 2),
      longtable = FALSE,
      center = "centering")


## Estimation of the proportions of semi-sincere and kick-out voters
## M = 2

propM2 <- dataM2 %>% 
summarise(ko_ub1 = mean(pref_runnerup == pref_topmplus1_third & vote_lastwin == 1),
					ko_lb1 = mean(pref_runnerup == pref_topmplus1_third & vote_lastwin == 1 & 
												vote_topmplus1_first == 0),
					ko_ub2 = mean(pref_runnerup != pref_topmplus1_third & vote_runnerup == 1),
					ko_lb2 = mean(pref_runnerup != pref_topmplus1_third & vote_runnerup == 1 & 
												vote_topmplus1_first == 0),
					ss_ub = mean(vote_topmplus1_first == 1),
					koss1 = mean(pref_runnerup == pref_topmplus1_third & vote_lastwin == 1 & 
												vote_topmplus1_first == 1),
					koss2 = mean(pref_runnerup != pref_topmplus1_third & vote_runnerup == 1 & 
												vote_topmplus1_first == 1)) %>%
mutate(ko_ub = ko_ub1 + ko_ub2,
			 ko_lb = ko_lb1 + ko_lb2,
			 koss = koss1 + koss2,
			 ss_lb = ss_ub - koss,
			 propD = 1 - koss - ko_lb - ss_lb,
			 propBD = 1 - ss_ub,
			 propCD = 1 - ko_ub,
			 prop_cross_pressured = ko_lb / (ko_lb + ss_lb))
propM2 <- format(round(propM2 * 100, 1), nsmall = 1)
propA <- paste0(propM2$koss, "\\%")
propB <- paste0(propM2$ko_lb, "\\%")
propC <- paste0(propM2$ss_lb, "\\%")
propD <- paste0(propM2$propD, "\\%")
propAB <- paste0(propM2$ko_ub, "\\%")
propAC <- paste0(propM2$ss_ub, "\\%")
propBD <- paste0(propM2$propBD, "\\%")
propCD <- paste0(propM2$propCD, "\\%")
propABCD <- paste0("100.0\\%")

## Table 1
rnameP_M2 <- c("", 
	             "",
	             "\\multirow{2}{*}{Kick-out Voting}",
						   "",
						   "Total")
resultP_M2 <- data.frame(rnameP_M2,
	                    c("", "", "Yes", "No", ""),
											NA,
											c("\\multicolumn{2}{c}{Semi-sincere Voting}", "Yes", propA, propC, propAC),
											NA,
											c("", "No", propB, propD, propBD),
											NA,
											c("\\multirow{2}{*}{Total}", "", propAB, propCD, propABCD))
captionP_M2 <- "Proportions of the Kick-out and Semi-sincere Voters ($M = 2$)"
latex(object = resultP_M2,
      title = "",
      file = "tables/experimentP_M2.tex",
      label = "tb_experimentP_M2",
      caption = captionP_M2,
      first.hline.double = FALSE,
      colheads = FALSE,
      rowname = NULL,
      col.just = c("l", rep("c", ncol(resultP_M2) - 1)),
      longtable = FALSE,
      center = "centering")

## M = 3
propM3 <- dataM3 %>% 
summarise(ko_ub1 = mean(pref_runnerup == pref_topmplus1_least & vote_lastwin == 1),
					ko_lb1 = mean(pref_runnerup == pref_topmplus1_least & vote_lastwin == 1 & 
												vote_topmplus1_first == 0),
					ko_ub2 = mean(pref_runnerup != pref_topmplus1_least & vote_runnerup == 1),
					ko_lb2 = mean(pref_runnerup != pref_topmplus1_least & vote_runnerup == 1 & 
												vote_topmplus1_first == 0),
					ss_ub = mean(vote_topmplus1_first == 1),
					koss1 = mean(pref_runnerup == pref_topmplus1_least & vote_lastwin == 1 & 
												vote_topmplus1_first == 1),
					koss2 = mean(pref_runnerup != pref_topmplus1_least & vote_runnerup == 1 & 
												vote_topmplus1_first == 1)) %>%
mutate(ko_ub = ko_ub1 + ko_ub2,
			 ko_lb = ko_lb1 + ko_lb2,
			 koss = koss1 + koss2,
			 ss_lb = ss_ub - koss,
			 propD = 1 - koss - ko_lb - ss_lb,
			 propBD = 1 - ss_ub,
			 propCD = 1 - ko_ub,
			 prop_cross_pressured = ko_lb / (ko_lb + ss_lb))
propM3 <- format(round(propM3 * 100, 1), nsmall = 1)
propA <- paste0(propM3$koss, "\\%")
propB <- paste0(propM3$ko_lb, "\\%")
propC <- paste0(propM3$ss_lb, "\\%")
propD <- paste0(propM3$propD, "\\%")
propAB <- paste0(propM3$ko_ub, "\\%")
propAC <- paste0(propM3$ss_ub, "\\%")
propBD <- paste0(propM3$propBD, "\\%")
propCD <- paste0(propM3$propCD, "\\%")
propABCD <- paste0("100.0\\%")

## Table D.1
rnameP_M3 <- c("", 
	             "",
	             "\\multirow{2}{*}{Kick-out Voting}",
						   "",
						   "Total")
resultP_M3 <- data.frame(rnameP_M3,
	                    c("", "", "Yes", "No", ""),
											NA,
											c("\\multicolumn{2}{c}{Semi-sincere Voting}", "Yes", propA, propC, propAC),
											NA,
											c("", "No", propB, propD, propBD),
											NA,
											c("\\multirow{2}{*}{Total}", "", propAB, propCD, propABCD))
captionP_M3 <- "Proportions of the Kick-out and Semi-sincere Voters ($M = 3$)"
latex(object = resultP_M3,
      title = "",
      file = "tables/experimentP_M3.tex",
      label = "tb_experimentP_M3",
      caption = captionP_M3,
      first.hline.double = FALSE,
      colheads = FALSE,
      rowname = NULL,
      col.just = c("l", rep("c", ncol(resultP_M3) - 1)),
      longtable = FALSE,
      center = "centering")
